#include <iostream>
#include <fstream>
#include "TDIE_cylinder.h"
using namespace std;
double I[M][N];
Space3D Rho[M];
double I_func(int k,double time)
{
	//if (time/dt == int(time/dt)) return I[k][int(time/dt)];
	//else return I[k][int(time/dt+0.5)]*(1-(int(time/dt+0.5)*dt-time)/dt)+I[k][int(time/dt-0.5)]*(1-(time-int(time/dt-0.5)*dt)/dt);
	//return I[k][int(time/dt+0.5)];
	return (I[k][int(time/dt+0.5)]+I[k][int(time/dt-0.5)])/2;
}
void Geometry()
{
	for (int m =0; m< M; m++)
		Rho[m] = Space3D(R * cos(Angle_distance * m), R * sin(Angle_distance * m) ,0);
}
double E_incident(int m,int n)
{
	double Gamma = 0;
	Gamma = 4/T * (C*n*dt - vt0 - Rho[m] * a_k);
	return E0*4 / (T * sqrt(pi)) * exp(-Gamma * Gamma);
}

double E_integral(int m,int n)
{
 	double temp = 0;
	for (int i =0; i< n; i++)
		temp += (E_incident(m,i) + E_incident(m,i+1)) * dt / 2;
	return temp;
}
double R_func(int m, int k, int l)
{
	Space3D temp = Rho[m] - Rho[k];
	return sqrt(temp.x*temp.x+temp.y*temp.y+temp.z*temp.z+(l*Delta_length)*(l*Delta_length));
}
double Alpha_func(int m, int k, int l)
{
	return Delta_length*Delta_length/R_func(m,k,l);
}
void Output()
{
	ofstream cur;
	cur.open("cur.txt");
	//for (int n=0; n<N; n++)
	//	cur<<E_incident(0,n)<<endl;
	for (int n=0; n<N; n++)
			cur<<I[0][n]<<endl;
}
int main()
{
	Geometry();
	for (int m=0;m< M; m++)
		I[m][0] = 0;
	for (int n=1;n< N; n++)
		for (int m=0; m< M; m++)
		{
			double temp = 0;
			for (int k=0; k<M ; k++)
			{
				if (k!=m)
				{
					if (n*dt - R_func(m,k,0)/C >=0)
						temp = temp + Mu/(4*pi) * I_func(k,n*dt - R_func(m,k,0)/C) * Alpha_func(m,k,0);
				}
				for (int l=1;;l++)
				{
					if (n*dt - R_func(m,k,l)/C >=0)
						temp = temp + Mu/(2*pi) * I_func(k,n*dt - R_func(m,k,l)/C) * Alpha_func(m,k,l);
					else break;
				}
			}
			temp = E_integral(m,n) - temp;
			I[m][n] = temp * pi / (Mu * log(1+sqrt(2.0)) * Delta_length);
			if (n>1) I[m][n-1] = 0.25*(I[m][n-2]+2*I[m][n-1]+I[m][n]);
		}
	Output();
	return 0;
}
